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One widely used technique for the construction of equihbrium models of stellar disks is 
based on the Jeans equations and the moments of velocity distribution functions derived 
using these equations. Stellar disks constructed using this technique are shown to be "not 
entirely" in equilibrium. Our attempt to abandon the epicyclic approximation and the 
approximation of infinite isothermal layers, which are commonly adopted in this tech- 
nique, failed to improve the situation substantially. We conclude that the main drawback 
of techniques based on the Jeans equations is that the system of equations employed is 
not closed, and therefore requires adopting an essentially ad hoc additional closure condi- 
tion. A new iterative approach to constructing equilibrium iV-body models with a given 
density distribution is proposed. The main idea behind this approach is that a model is 
first constructed using some approximation method, and is then allowed to adjust to an 
equilibrium state with the specified density and the required parameters of the velocity 
distribution remaining fixed in the process. This iterative approach was used to construct 
isotropic, spherically symmetric models and models of stellar disks embedded in an ex- 
ternal potential. The numerical models constructed prove to be close to equilibrium. It is 
shown that the commonly adopted assumption that the profile of the radial velocity dis- 
persion is exponential may be wrong. The technique proposed can be applied to a wide 
range of problems involving the construction of models of stellar systems with various 
geometries. 
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1 Introduction 



In studies of the dynamic evolution of galaxies using the results of numerical A^-body 
simulations, it is very important to correctly specify the initial equilibrium state of the 
stellar system. Two different approaches are employed to achieve this. The first approach 
uses the kinetic equation (coUisionless Boltzmann equation) and the second approach 
deals with the equations of stellar hydrodynamics (Jeans equations). Both approaches 
have their advantages and disadvantages. 

The kinetic approach is based on the use of the phase-space distribution function 
(DF) and the Jeans theorem, which states that any function of integrals of the motion is 
a solution of the stationary coUisionless Boltzmann equation [1]. For spherically symmet- 
ric stellar systems with isotropic velocity distributions, any function of the form f{E), 
where E is the specific energy, describes an equilibrium gravitating system. If spherically 
symmetric models with a density distribution profile p{r) that corresponds to the obser- 
vational data and a potential $(r) derived from the Poisson equation are applied to real 
objects (elliptical galaxies and various subsystems of spiral galaxies, such as the bulge and 
halo), determining the distribution function f{E) reduces to solving an integral equation 
derived via the Abel transform (Eddington's formula; see. e.g., [I]). Analytical solutions 
of this equation are known only for special classes of modeldj. Researchers often prefer 
a different approach, proceeding from a particular form of f{E) and integrating it over 
velocity space to obtain the density distribution and then the potential. One can then 
choose from among the broad class of analytical models constructed in this way those 
whose density profiles are closest to the density profile of the real system. 

The transition to anisotropic or axisymmetric models, e.g., stellar-disk models, drasti- 
cally complicates the problem. Describing a disk system requires the use of a distribution 
function of the form f{E,Lz), where Lz is the angular momentum of a particle with 
respect to the symmetry axis z. A number of models are known to have analytical for- 
mulas simultaneously for f{E,Lz), p, and $ (simple models include the Kalnajs [5J disk, 
however, this is not the only one). Such models can be used to study the properties of 
flat systems only as a flrst approximation. These are usually 2D models and their radial 
proflles usually differ strongly from the exponential form typical of spiral galaxy disks. 

As for 3D axisymmetric models, constructing such models using the distribution func- 
tion f{E, Lz) also faces another serious barrier. In systems represented by such distri- 
bution functions, the radial velocity dispersion should be equal to the velocity dispersion 
parallel to the rotation axis (see, e.g., [1]), which is inconsistent with observational data - 
at least with data for our Galaxy based on measurements of stellar velocity dispersion in 
the solar neighborhood (see, e.g., [6]). In the kinetic approach, axisymmetric models with 
anisotropic radial and vertical velocity distributions can be represented by a function of 
the form f{E, Lz, I3), where I3 is the third integral of the motion. The general expression 
for I3 is unknown. In sufficiently cool stellar disks, where the stellar velocity dispersion is 
small compared to the velocity of the disk's rotation about the symmetry axis, the energy 
of vertical oscillations Ez = $(-R, z) — $(-R, 0) + |f ^ is a well-conserved quantity (along 
almost circular orbits). This can be used as a third integral of motion when constructing 

"'^For example, the distribution function corresponding to the popular dark halo model - the Navarro- 
Frenk- White 2] model - can be determined only numerically [3], [1]. 
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the DF for thin stellar disks with density distributions pdisk{R,z), close to the observed 
exponential law [7], [H]. Kuijken and Dubinski [7] and Widrow and Dubinski [S] also 
describe a procedure for correcting the DF in the case of multicomponent galaxy models. 
These authors actually solve a more general problem - the construction of the distribu- 
tion function for a disk-galaxy model consisting of several self-consistently components: 
an exponential disk, bulge, dark halo, etc. However, even more specific problems, such as 
a construction of an equilibrium model for a disk with a realistic density profile, which is 
embedded in the external potential, is not easy to produce as an initial configuration for 
A^-body simulations. 

The second approach, which is based on computing the moments of the equilibrium 
particle velocity distribution function, uses the Jeans equations. In the case of construct- 
ing equilibrium stellar-disk models with realistic density profiles, this approach usually 
consists of various modifications of the technique described by Hernquist [U]. 

This technique has the advantage that it is relatively simple and makes it possible to 
construct a model that is more or less close to equilibrium for a given density distribution 
profile pdisk{R,z) and external potential $ext(-R, -z), under reasonable assumptions about 
the velocity distribution function. However, this technique has an important drawback: 
the resulting A^-body model is often far from equilibrium. This primarily concerns mod- 
els with small masses for the spheroidal components (bulge and halo). We analyze the 
technique of Hernquist in more detail in the next section. 

In this paper, we offer a new iterative method for constructing equilibrium models of 
stellar disks embedded in an external (rigid) potential, which we compare to the widely 
used method based on the moments of the distribution function. Our iterative models are 
very close to equilibrium and have the specified density profile. The iterative approach 
has a broader scope than the particular problem we are solving, and can be used to both 
construct models with a different geometry (e.g., spherically symmetric models with a 
given mass distribution and anisotropy of the random motions) and construct equilibrium 
models of selfconsistent multicomponent stellar systems. 

2 Approach based on the moments of the distribu- 
tion function 

If a stellar disk with density Pdisk(-R, ^) and an external potential $ext(-R, ^) produced, 
e.g., by the halo and bulge, are axisymmetric, the Jeans equations for the first and second 
moments of the velocity distribution function for the disk particles can be written in the 
form P, [in] 
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Pdisk dR 
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where is the mean azimuthal velocit}@; ctr, a^, and cr^ are the velocity dispersions 
in the radial, azimuthal, and vertical directions, respectivel}{§; $tot = '^'ext + '^'disk is the 

total potential produced by all components of the system: and Vc = R ^ is the circular 

oR 

velocity. Here, we have omitted the dependences of parameters on the coordinates R and 
z in the cylindrical coordinate system for simplicity. 

Equations ([T]) assume the absence of regular motions in the radial and vertical direc- 
tions. We also assume that the axes of the velocity ellipsoid are directed along the axes of 
the cylindrical coordinate system; this means that second moments of the form vrVz are 
equal to zero. This is a reasonable assumption for the galactic plane (based on symmetry 
considerations). However, outside the galactic plane, the velocity dispersion ellipsoid is 
inclined so that the equality vrVz = breaks down, although substantial deviations 
appear only in the central regions of the disk. 

The system ([T]) consists of three equations and has four unknowns: v^, an, a^p, and 
cr^. To solve such a system, one should make some additional assumption, and this is the 
major drawback of methods of constructing equilibrium models for stellar disks based on 
the Jeans equations. 

One of the most popular techniques for constructing equilibrium (or, more precisely, 
close-to-equilibrium) A^-body stellar-disk models based on the Jeans equations is described 
by Hernquists p. Many authors (see, e.g., [TT], [12], [13], [E], [E]) have applied this 
technique with various minor modifications. It is usually applied to three-dimensional 
disks with exponential density profiles: 

'^^^^ exp ( — ] sech^ f ~ ) , R ^ Rn 



PdM z) = { Anh^zo ^ V hj \zoJ ' " '"'^ ' (2) 

, R> -Rmax . 



Here, h is the exponential disk scale, Zq is the vertical scale length, Mdisk is the total 
mass of the dislo, and -Rmax the disk radius (truncation radius). This density profile 
approximates well the observed profiles of real spiral galaxies [16j. 

Further, the conditions under which the Jeans equations ([1]) were derived should be 
supplemented with the following additional assumptions [0]. 

1. All four moments {v^, aR, a^p, and az) are independent of z and depend only on the 
cylindrical radius R; i.e., the disk is isothermal in the z direction. 

2. The epicyclic approximation is valid (random velocities of stars are small compared 
to the velocity of rotation). In this case, the mean azimuthal velocity in the second 
equation of ([T]) can be replaced by the circular velocity (i.e., we can substitute Vc 
for v^). 

3. The last equation of ([1]) is often rewritten in the approximation of infinitesimal 
isothermal layers [T7]. In this case, the contribution of the external potential to 



^An over line denotes averaging. 

^In accordance with the usual astronomical usage, the dispersion denotes the standard deviations of 
the distribution function. 

^Formally, Mdisk is the total mass of the disk at i?inax = oo. 
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the total potential $tot is often neglected. This yields the relation o"^ = vrGSdisk^^o, 
where Sdisk is the surface density of the disk. 

4. The velocity distribution function has the Schwartzschild form, i.e., the velocity 
components along each of the three axes of the cylindrical coordinate have normal 
distributions. 

These are standard assumptions, are believed to be appropriate for real galaxies p. 

The system ([T]) is closed under the additional assumptions: cr^ oc exp{—R/h). It 
is convenient to specify the coefficient of proportionality via the Toomre parameter Qt 
[T8] at some radius -Rref- This parameter characterizes the degree of disk heating, or 
the margin of stability against the growth of perturbations in the disk plane. For a 
stellar disk with an exponential density profile in R, the relation cr^ oc exp{—R/h) means 
that o"|. is proportional to the surface density Sjisk- Together with the isothermal layer 
approximation, this also yields a dependence of the form oc (7^. Observational data for 
our Galaxy [19], [20] are consistent with the dependence cr^ oc exp{—R/h). It is believed 
from general considerations that this relations is also appropriate to other spiral galaxies 
[21] . We examine the validity of this assumption below. 

As a result, one can construct for a given p^isk and $ext a one parametric family of 
models parametrized by a quantity that characterizes the degree of disk heating, or the 
fraction of the kinetic energy of the disk contained in random motions. 

The deficiency of this approach is that the constructed stellar disks prove to be not en- 
tirely in equilibrium. Although the models rapidly adjust to equilibrium, this adjustment 
may complicate analyses of the results of A^-body simulations, such as those used to study 
instabilities of the stellar disk. In the process of this adjustment, a characteristic ringlike 
density wave forms and propagates from the center, as is illustrated by the results of our 
numerical simulation^ in Fig. [1] (see [H], [15] for a more detailed description of the tech- 
nique used for the numerical simulations). Kuijken and Dubinski [7] also pointed out a 
similar effect. The most nonequilibrium models are hot models (with large QT(-Rref) ~ 2) 
and models without halos or with low mass halos. The only nearly equilibrium models are 
those with fairly massive halos (such as the models described by Hernquist 0, where the 
mass of the halo within four exponential disk scales was greater than five disk masses). 

The nonequilibrium nature of the models constructed using the Hernquist technique is 
due to the underlying assumptions of this technique. Our attempt to refine this approach 
by abandoning the epicyclic and isothermal layer approximations did not result in any 
substantial improvement of the models. 

For example, the solution of the equation of vertical equilibrium (the third equation 
of ([T])) and of the equation derived from it in the isothermal layer approximation yield 
similar az{R,z) in many cases, and have little effect on the extent to which the models 
are nonequilibrium. The results of the ensuing computations illustrate this conclusion. 
Since the equation of vertical equilibrium contains only one unknown, cr^, when p{R,z) 

^In all numerical simulations that used stellar-disk model (2), we set G ^ 1, h = 3.5, and Mdisk = 1- If 
needed, the results of our simulations can be interpreted in a dimensional system of units (one of several 
possible) in which the unit of length is Ru = 1 kpc and the unit of mass is M„ = 8 • IO^^Mq. The units 
of time and velocity are then tu ~ 1.67 Myr and « 587 km/s. 
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Figure 1: Adjustment to equilibrium of the model constructed using the Hernquist tech- 
nique (see Section [2]). The top snapshots show a face-on view for several times, with the 
shades of gray corresponding to the logarithm of the particle number per pixel. The plots 
at the bottom show the radial distribution of the number of particles (counts made in con- 
centric cylindrical layers). A characteristic wave propagating from the center can be seen. 
The exponential disk model ([2]) is used with parameters h = 3.5, zq = 1, Mdisk = 1, and 
-Rmax = 14. The external potential has the form of a Plummer sphere (HI) with a^i = 15 
and Mpi = 4. For these parameters, the total fractional mass of the spherical component 
within four exponential disk scales is about 1.3 (i.e, Msph(4/;,)/Mciisk(4/i) ~ 1.3). The 
Toomre parameter is QT(-Rref) = 1-5, where i?ref = 8.5. 150 time units approximately 
corresponds to one disk rotation period aX R = 8.5. The number of particles, smoothing 
parameter, and integration step are = 25000, e = 0.02, and dt = 0.01, respectively. 
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is specified, this equation can be solved independently of the other equation^. As is 
evident from Fig. [21 the approximation of infinite isothermal layers breaks down only if 
the external potential in the system is represented by a very massive dark halo (the "disk 
+ heavy halo" curves) or varies strongly over a vertical disk scale height, e.g., when the 
system has a bulge (the "disk + bulge" curves). If the system has no external potential 
at all (the "disk" curves) or is represented by a small halo with a mass comparable to the 
mass of the disk (the "disk + light halo" curves), then the two methods yield virtually the 
same dispersion (except for the centermost regions, where the approximation of infinite 
isothermal layers obviously breaks down, although, even there, the differences are not 
very large). 

We believe that the nonequilibrium state of the models is due primarily to the adopted 
additional assumption concerning the dependence of crji on R. 

Note that a number of authors have adopted other assumptions, different from the 

proportionality (t|. oc exp{—R/h), or equivalently ctr oc cTx, as the missing equation in the 

system of Jeans equations. For example, Athanassoula and Misiriotis [13] adopted the 

additional assumption that the Toomre parameter Qt is independent of radius. Revaz 

u 

and Pfenniger [IT] assumed that an oc cr^ — , where k is the epicyclic frequency and u is 

the frequency of the vertical oscillations, (z/^ = ° ). They computed the coefficient 

of proportionality as in the original technique of Hernquist, by specifying the Toomre 
parameter at some radius. 

So far, we have no grounds to prefer any particular additional condition to supplement 
the Jeans equations. However, by taking into account all the moments up to the sixth 
inclusive, and if that all the moments beginning with the second are small compared to 
the circular velocity, this problem can be solved without invoking additional assumptions 
[23], [21] • Unfortunately, due to its complexity and awkwardness, this solution has never 
been used to construct A^-body models of stellar disks. 



3 Dependence of the radial velocity dispersion on the 
cylindrical radius (dependence of on R) 

The assumption that aR oc az is now generally accepted, or, in any case, it is used to 
analyze and interpret the observation data [21], [25], [26]. However, is this assumption 
justified? 

Thus far, the observed dependence of aji on R has been obtained only for the Milky 
Way. Based on a derivation of this dependence from observations of K giants, Lewis and 
Freeman [19] concluded that aj^ oc exp{—R/h). In the isothermal layer approximation, 
this means that cr/j oc o"^. However, we must emphasize an important point. The errors in 

^For example, Revaz and Pfenniger [TT], Khoperskov et al. [H], and Bahcall used the solution of 
this equation in analyses of the equilibrium in the vertical direction and the construction of equilibrium 
models for stellar disks. Note that this equation can be solved in two ways. As we did here, one can 
specify the density profile and find az{R,z). Alternatively, iJz{R,z) can be specified - for example, by 
assuming that cr^ is independent of z ^22, - in order to determine the vertical density profile. Khoperskov 
et al. [12j implemented the latter approach. 
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Figure 2: Dispersion of the vertical stellar velocities, az, computed both in the approx- 
imation of infinite isothermal layers and without this approximation for various external 
potentials. The top left plot shows the computed dispersion in the disk plane {z = 0), 
the top right plot the dependence of the dispersion on 2; at i? = 0, and the bottom right 
plot the dependence of the dispersion on z at R = 3.5. In this case, we used a model of a 
relatively thin exponential disk ([2]) with h = 3.5, zq = 0.3, Mdisk = 1, and -Rmax = 00. The 
"old" curve shows the dispersion computed in the approximation of infinite isothermal 
layers. All other curves were constructed without this assumption and the corresponding 
models differ in the adopted external potentials: the "disk" model was computed with no 
external potential, the "disk + light halo" model with an external potential in the form 
of Plummer sphere (HI) with Opi = 3.5 and Mpi = 1, the "disk + heavy halo" curve with 
an external potential in the form of the Plummer sphere with Opi = 3.5 and Mpi = 5, and 
the "disk + bulge" curve with an external potential in the form of a Hernquist sphere 
with Ohr = 0.5 and Mhr = 0.1. 
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the dispersion obtained by Lewis and Freeman [19] are fairly large (on the order of 10%), 
as are the errors in the stellar distances, which were derived indirectly. Accurate (errors of 
the order of one percent) measurements of are available only in the solar neighborhood, 
based on HIPPARCOS data |6j. Therefore, for the Galaxy, we can only conclude that 
aR is approximately proportional to az- It seems possible that, if the dependence of a^p 
on R for the Galaxy could be measured with the same accuracy as we now know the 
dependence for ct/j, it would also fit the relation oc exp{—R/2h) oc a^. This situation 
should improve when new astrometric satellites (in particular, GAIA) are launched. At 
least for the Galaxy, the dependence of the moments of the velocity distribution function 
on R will then be known more accurately. 

It is not possible to obtain aR{R) for other galaxies directly from observations, which 
can yield only the hne-of-sight velocity, fios, and line-of-sight velocity dispersion, aios- For 
edge-on galaxies, aios depends on both and a^, whereas, for intermediate- inclination 
galaxies, all three components of random velocities {a^, a^, a^) contribute to aios- How- 
ever, we can use the Jeans equations ^ to derive the unknown dependence on R 
from the observed quantities fios and aios- Gerssen et al. [27], [26], Shapiro et al. [28] . 
and Wesfallet al. [29] used such considerations to infer the velocity dispersions in three 
directions from observational data for inclined galaxies. The drawback of the technique 
they used is that it involves an a priori assumption about the form of the dependence of 
aR and on R: 



where a, a^fi, and aR^, are parameters determined during the reduction of the data. 
These authors pointed out that there are no grounds to believe that the scale parameter 
a should be the same for and aR, but the quality of the observational data prevents 
independent determination of these two scale parameters. 

It follows that the quality of the observational data available so far prevents the deter- 
mination of the dependence of aR on R for external galaxies, and available observations 
can only yield estimates of the velocity dispersion for the entire galaxy. 

We also note the two semi-theoretical papers [3^ and [31], whose authors used inde- 
pendent methods to conclude that aR oc is not valid for our Galaxy. 

4 Iterative method for constructing equilibrium 
models 

A^-body simulations sometimes use the following method to specify the initial conditions. 
The initial conditions are specified (in some way) to be close to equilibrium. The system 
is then given some time to adjust to a new equilibrium state, which is used as the initial 
state for the A^-body simulations (Sellwood and Athanassoula [32] and Barnes [33] used 
this approach with some modifications). The drawback of this technique is that it is 
difficult to build a close-to-equilibrium model with a given density profile. 

We developed a method that is essentially a logical extension of the above technique. 
The main idea of our method is to allow the system to adjust to the equilibrium state 



az{R) 

MR) 




a^fi exp(-i?/a) , 



(3) 
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while fixing the density profile. We achieve this via the following general algorithm for 
the iterative method. 



1. Use some approximate method to construct a close-to-equilibrium iV-body model 
with agiven density distribution (i.e., with a given particle distribution function in 



2. Allow the model to evolve for a short time. 

3. Construct a model with the same velocity distribution as that of the evolved model, 
but with the required density profile (the initial density profile). Note that, if we 
have certain constraints on the particle velocity distribution (e.g., if we want to 
construct a Plummer sphere model with an isotropic velocity distribution), we must 
correct the particle velocity distribution to satisfy these constraints. 

4. Return to item 2. Stop iterations when the velocity distribution ceases to change. 

As a result, we obtain a close-to-equilibrium A^-body model with the given density profile. 

Practical implementation of the iterative method is somewhat more difficult. The 
main problem arises at the third stage - the construction of a model with the same 
velocity distribution as for the slightly evolved model of the previous iteration step. In 
the ideal case, we would have to obtain the particle velocity distribution at each point of 
the system. This is, naturally, impossible, given the available number of particles in the 
A^-body simulation. However, this problem can be resolved if we make some simplifying 
assumptions. 

The next section describes how an iterative method can be used to construct a spheri- 
cally symmetric, equilibrium, isotropic model with a given distribution function (a Plum- 
mer sphere). In the problem at hand, this enables us to verify that the iterative method 
indeed permits the construction of an equilibrium model. This realization of the tech- 
nique can be used to construct spherically symmetric, isotropic-velocity models with other 
density profiles and, with small modifications, even models with anisotropic velocity dis- 
tributions. Below we apply the iterative method to construct equilibrium A^-body models 
for stellar disks. The realization of the third part of the algorithm in the case of spheri- 
cally symmetric models differs from that for disk models, but the main idea remains the 
same. 

4.1 Equilibrium, isotropic, spherically symmetric models 

To test the idea of the iterative method, we will use it to construct an equilibrium model 
for a case when the equilibrium distribution function is known a priori: an equilibrium, 
isotropic model of a Plummer sphere. The isotropy means that there is no special direction 
in velocity space, i.e., the velocity distribution function depends only on the speed. The 

'''This initial model does not have to be close to equilibrium. In the next section, we use a cold model 
(with zero velocities) as an initial approximation to construct an equilibrium Plummer sphere. 
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equilibrium distribution function for this model is known (see, e.g., [T], p. 223), and the 
potential for this model has the form 




(4) 




2 ' 



where Mpi is the total mass for the model and Opi is the scale length (the Plummer model 
has about 35% of its mass inside the radius Opi). 

We now implement the third part of the algorithm of the iterative method (the "trans- 
fer" of the velocity distribution function) as follows. We take our slightly evolved model, 
from which we plan to copy the velocity distribution function. We subdivide this model 
into spherical layers containing approximately the same number of particles and construct 
the distribution of the particle speeds v in each of these spherical layers. 

We do this in the usual way, by determining the range of variation of v and subdividing 
this interval into some number of bins. We then compute the number of particles falling 
in each bin. In the limit of an infinite number of particles and infinite number of bins, the 
resulting histogram gives the distribution function (not normalized to unity). We assumed 
that the number of particles in a bin is equal to the value of distribution function at the 
center of the bin. We can then approximate these values by some function and use this 
as the distribution function. We used several types of approximation: piecewise-linear, 
cubic spline, and least-squares fits using functions of the form exp(P(x)), where P{x) is a 
polynomial. The type of fit had virtually no effect on the result of the iterative algorithm. 

This yields the distribution function of velocity modulus in each of the spherical layers. 
We then constructed the model for the next iteration step as follows. We specified the 
positions of the particles in accordance with the given density profile. We then determined 
the spherical layer to which each particle belongs, and used the "selection-rejection" 
method with the distribution function for the given spherical layer to determine the speed, 
choosing the direction for the velocity as random. 

The velocity distribution function for the new model is isotropic and spherically sym- 
metric (here, we mean that the velocity distribution function is symmetric in coordinate 
spac^. We did not simply transfer the velocity distribution function, but adjusted it to 
make it strictly isotropic and spherically symmetric. 

We now apply the algorithm described above to construct an equilibrium model for the 
Plummer sphere. We use virial units (G = 1, the total mass of the model is Mpi = 1, and 
the total energy is -E = — 1/4). In iterative models, we truncated the density distribution 
at Tmax = 5 (the sphere of this radius contains about 98% of the mass of the modeled 
system) . 

Let us consider two models constructed using iterative method: 

• model "li" with a piecewise-linear approximation for the distribution function; 

• model "epol," in which the distribution function is approximated by functions of the 
form exp(P(x)), where P{x) is a polynomial (we used a sixth-degree polynomial). 

^Isotropy can be considered to be spherical symmetry in velocity space. 
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We constructed both models in 120 iterations, with the first 100 having relatively low 
accuracy. The number of particles is = 10^, the number of spherical layers used to 
subdivide the system rir = 50, and the number of bins used to construct the distribution 
function of the speed v in each spherical layers, = 21. We performed the last 20 
iterations with relatively high accuracy: = 5 ■ 10^, rir = 200, = 31. The time 
step for each iteration is = 1 (which approximately corresponds to the crossing time of 
the system's core in the adopted units). Our starting point was a cold model with zero 
velocities. 

Note that the initial model in our computations is far from equilibrium; however, the 
iterations converge to close-to-equilibrium models. Figure [3] demonstrates the convergence 
of the iterations during the construction of model li. It is clear that the iterations converge 
to the equilibrium model. Departures from equilibrium are evident only at the periphery 
of the system. For example, the velocity dispersion, a^, for model li differs appreciably 
from the value given by the equilibrium model at r > 3 (Fig. [H]). However, the system 
has only about 5% of its mass located outside the sphere of radius r = 3. 

Let us now experimentally test the equilibrium of our models by comparing them 
with two other models. The first one is constructed using the equilibrium distribution 
function as an example of the closest-to-equilibrium model (model "DF"), and the second 
is constructed using another approximation method, based on the following main idea. 
We can compute the mean speeds, v, and velocity dispersion, a^, at each radius of the 
equilibrium model for the Plummer sphere: 

V ^ 0.665-^-<l>pi(r) , 

(5) 

a, ^ 0.240y/-<l>pi(r). 

The numerical coefficients in are determined in virial units by integrating the equilib- 
rium distribution function (in virial units, Mpi = 1 and a^i = 37r/16). 

An approximate model can be constructed assuming that the speed distribution func- 
tion is Gaussian with the moments ([5]); we will call this as the "gauss" model. 

Let us compare the initial stages of the evolution of the four A^-body models for the 
Plummer sphere (li, epol, DF, and gauss) for A^ = 50 000. We chose the integration step 
and smoothing parameter for the potential e in A^-body simulations in accordance with 
the recommendations derived in [M] {dt = 0.002, e = 0.007). 

We analyzed the results using the and A^, values introduced and determined as 
follows (we used similar quantities in our earlier paper [31]). 

The quantity A^ characterizes the spatial deviation of the particles at time t from 
the initial distribution (at time t = 0). We computed this by subdividing the model into 
spherical layers, computing in each spherical layer the difference between the numbers 
of particles at time t and at the initial time, and then computing A^ as the sum of the 
absolute values of these differences, normalized to the total number A^ of particles in the 
system. The thickness of the layer and maximum radius were 0.1 and 2, respectively. 

Note the very important point that the value of A^ for two random realizations of any 
model is not zero. This quantity has appreciable natural noise due to the finite number of 
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Figure 3: Convergence of iterations in the construction of model li. The upper plots show 
the dependence of the mean speed, v, on r, while the lower plots show the dependence of 
the dispersion of the speed, a^, on r for several iterations. The left-hand plots show these 
quantities for the 1st, 10th, 50th, 100th, and 120th iterations, while the right-hand plots 
correspond to the last five iterations. The solid bold curve shows v and for equilibrium 
model (these quantities can be computed by integrating the equilibrium distribution 
function) . 
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particles. We estimated this noise level by computing the average value of for a large 
number of pairs (200) of random realizations of the model at the initial time. 

We computed A^, in the same way, but for the distribution of particles in velocity 
space. For the results reported here, the width of the layer and maximum velocity were 
0.1 and 1.5, respectively. 

Figure H] shows the time dependences of A^ and A^, for our four models. It is obvious 
that the values of A,, and A^ for models DF, li, and epol are close to their natural noise 
levels. At the same time, the A^ value for the gauss model appreciably exceeds the 
noise level; A„ also appreciably exceeds the noise level, although we note that the gauss 
model adjusts to the equilibrium state in about 15-20 time units (which approximately 
correspond to the crossing time for the entire system), and this regime has virtually the 
same profile in the equilibrium state and at the initial time. 

We can draw the following general conclusion. The models constructed using the itera- 
tive method (li and epol) are essentially in equilibrium. They conserve their parameters as 
well as the model constructed using an equilibrium distribution function (model DF). The 
iterative models behave much better than the model constructed in the approximation of 
a Gaussian velocity distribution (model gauss). 

Our models exhibit deviations from equilibrium at the very edge of the system, at r > 3 
(as we would expect; Fig. [S]). On the whole, these models are very close to equilibrium. 
We emphasize again that we used cold models, which are very far from equilibrium, as 
the initial models for our iterations. However, through the iterations, these models came 
to virtually equilibrium states. Thus, the iterative approach has fully validated itself. In 
the next section, we apply this approach to more complex systems. 

4.2 Equilibrium models of stellar disks 

4.2.1 Implementation of the iterative method to construct an equilibrium 
model for a stellar disk 

Let us apply the iterative method described above to construct an equilibrium (A^-body) 
model for a stellar disk with a density distribution pdisk{R,z) embedded in an external 
potential $cxt(-R, -z)- We would expect a family of equilibrium disks to exist for a given 
Pdisk and $ext- This must be at least a one-parameter family characterized by the fraction 
of kinetic energy contained in random motions. This is where the construction of an 
equilibrium stellar disk differs from the construction of an equilibrium isotropic Plummer 
sphere - we must construct an entire family and not just one model. 

The main difficulty in the iterative method is how to transfer the velocity distribution 
function (step 3 - see the beginning of Section Hj). This is done using the moments of the 
distribution function and certain assumptions about the form of the velocity distribution 
function (here, we have some similarity to the Hernquist method). We did this as follows. 

We took the disk model from which we planned to "copy" the velocity distribution 
function (a slightly evolved model from the previous iteration step). We then subdivided 
the model into concentric cylindrical layers of infinite length along the z axis. We com- 
puted in each layer the four moments of the distribution function v^, a^, a^, and a^, and 
used these moments to specify the velocities in the model for the next iteration step. We 
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Figure 4: Time dependences of and A^, for four A'^-body models of the Plummer sphere. 
The horizontal line shows the noise level for A^ and A„. The initial distribution for the 
gauss model is far from equilibrium (e.g., from the velocity distribution in the DF model). 
We therefore computed A^ for the gauss model in two different ways: first, as for all the 
other models, by calculating the difference between the velocity distributions at the initial 
time and at time t (the curve "gauss"), and second, by calculating the difference between 
the velocity distribution at time t and the equilibrium velocity distribution, using the 
initial velocity distribution for the DF model (the curve "gaussl"). 
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assumed that the velocity distribution has a Schwartzschild form except for one feature: 
we removed all particles capable of escaping beyond the disk boundary. 

Let us explain this last condition. If we have specified an axisymmetric potential $tot 
(in our case $tot = '^'disk + '^'ext) and a particle is specified to have the coordinates R, z 
and velocities vr, f<^, v^-, this particle can reach the radius i?max (the disk boundary) if 
and only if the following inequality [l](p. 117) is satisfied: 

0.5ivl + Vl + VD + $tot(i?, Z) > $tot(i?max, 0) + . (6) 

We then searched for the parameters of this so "truncated" Schwartzschild distribution 
{v'^, cr^, (T^, a'^), that make its moments equal to the specified values (-u^, ur, a^, cx^). 
In some cases, this is impossible in the peripheral regions of the galaxy. We then used 
in these regions an "energy-truncated" Schwartzschild distribution instead of the radius- 
truncated distribution, i.e., we used the Schwartzschild distribution without particles with 
positive total energies (particles capable of escaping to infinity). 

This yields a velocity distribution function in each cylindrical layer, thereby solving the 
problem of transferring the velocity distribution function. Another technical point is the 
following. In each cylindrical layer, we have a distribution function with the parameters 
v'^, a'j^, a'^, o'^. We referenced the value of each of these parameters to the value at the 
midpoint of the cylindrical layer and interpolated them using piecewise-linear functions 
to obtain a velocity distribution function that is continuous in radius. 

We adopted the following important simplifying assumptions during the transfer of 
the velocity distribution function (that may fail in real stellar disks of spiral galaxies): 

• the disk is isothermal in the z direction; 

• by adopting the Schwartzshild velocity distribution, we implicitly assume that 

= 0. 

We wish to emphasize two important points. First, the particular implementation of 
the iterative method used does not matter - it is the idea behind the method that is 
important. Second, the main test of any method for constructing an equilibrium A^-body 
model should, naturally, be experimentally verifying the extent to which the model is 
close to equilibrium. 

We used the following model versions as initial models for the iterations. 

• Models constructed using the Hernquist technique (see Section [2]). 

• Cold models in which the particles move in strictly circular orbits. This model is 
theoretically in equilibrium, but is highly unstable, even on time scales of the order 
of the iteration time. Therefore, this model becomes slightly heated in the iteration 
process. 

• A "shifted" model. Let us assume that we have an equilibrium disk model con- 
structed using the iterative method. In our technique, this model is determined 
by the four moments of the velocity distribution function in each cylindrical layer 
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{v^, aR, a^, a^). We construct a new model with -u^ = cu^, a'^ = ^ctr, cr'^ = -^cr^, 
a'^ = a^, where c is the shift coefficient, which is close to unity, the primed moments 
refer to the shifted model, and the unprimed moments refer to the old, unshifted 
model. 

• Other versions. For example, a model constructed for a different external potential 
or for a disk with other parameters. 

The iterations for all these versions converged and yielded various models. However, 
as expected, all these models form a one-parameter family (Figs. [5] and El see below for a 
more detailed description of these models). 

Before starting our discussion of the properties of the iterative models, let us point 
out a few technical details. We must usually construct either several models in a family 
with the same pdisk and $cxt or some particular model of this family. One example might 
be a model with a given kinetic energy fraction contained in the random motions of stars 
(in this case, we can fix the Toomre parameter at some radius). The convergence of the 
iterations is faster the closer the initial disk is to equilibrium. Therefore, to save CPU 
time, it is advisable to use the following "tricks." We can initially perform the iterations 
with lower accuracy (fewer particles and a coarser dividing of the model into cylindrical 
layers during the transfer of the particle velocity distribution), and then carry out several 
iterations with higher accuracy. If we have already constructed one equilibrium model 
from a family, all the other models of this family can conveniently be constructed using 
shifted models (described above, where we enumerated possible initial models for the 
iterations). Such shifted models turn out to be much closer to equilibrium than, for 
example, models constructed using the Hernquist technique, and the iterations converge 
appreciably faster. 

To accelerate the convergence of the iterations and construct models with a fixed 
fraction of their kinetic energy contained in random motions, we can fix this energy 
fraction at each iteration (i.e., fix the degree of heating of the model). Instead of fixing the 
amount of energy contained in random motions, we can fix the amount of energy contained 
in regular motions. We fixed this latter quantity via the total angular momentum with 
respect to the z axis: 

N 
i=l 

where rrii, v^j, and Ri are the mass, azimuthal velocity, and cylindrical radius of the zth 
particle, respectively. 

We fixed Lz and, at each iteration step, having constructed the new model (with 
the same velocity distribution as in the slightly evolved model of the previous iteration 
step), we adjusted the azimuthal velocities of the particles to make the total angular 
momentum of the system exactly equal to the given L^. We did this as follows. Let 
be the given angular momentum and L'^, the current angular momentum. We specified 
the new azimuthal velocities of the particles using the relations 

= V + ^jv ' (8) 
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Figure 5: First family of models constructed using the iterative method. Dependences 
on cylindrical radius R are shown for four moments of the velocity distribution [v^p, a^, 
a^p, and a^; two top rows), and the Toomre parameter Qt and ratio cTz/ctr (bottom row). 
The cross indicates the model whose equilibrium test is shown in Fig. [9l The various 
curves in the plots correspond to various models in this family. 



18 




19 



where v'^^ is the current azimuthal velocity of the ith particle, v^i the adjusted azimuthal 
velocity of this particle, and Wi a weight that determines how we add the angular mo- 
mentum to the particle. In our experiments, we set Wi = Riirii. 

Iterations with fixed Lz converged to models of the same family as the iterations 
without fixed L^. However, the iterations with fixed Lz converged much faster. Fixing Lz 
is helpful, since it results in the convergence of the iterations to a strictly defined model 
from a family of iterative models with given p^isk and $ext; irrespective of the initial model. 



4.2.2 Properties of the iterative disk models 

Let us consider as examples two families of models constructed using the iterative method. 

Figure shows the first family. We adopted disk model ^ with h = 3.5, zq = 1, 
^disk = 1; and -Rmax = 14, and modeled the external potential using the Plummer 
sphere @ with Opi = 15 and Mpi = 4. For these parameters, the fractional mass 
of the spherical subsystem within four exponential radii is equal to about 1.3 (i.e., 
Miph(4/i)/Mdisk(4/i) ~ 1.3). Note that the parameters of the disk and external poten- 
tial are the same as in the model constructed using the Hernquist technique (Fig. [T]), 
which was used as an example of a not-quite-equilibrium model. The time for a single 
iteration, ti = 15, corresponds to approximately one-tenth of a disk rotation at a distance 
of two exponential disk scales. The thickness of the cylindrical layers into which we sub- 
divided the model when transfering the velocity distribution function is (Ir = 1. We used 
N = 10^ particles. 

Figure [U] shows the second family. It was produced using an initially thinner disk 
model with h = 3.5, Zq = 0.3, M^isk = 1, and i?max = 14, modeling the external 
potential using the Plummer sphere (jlj) with Opi = 3.5, Mpi = 0.88 and the potential of a 
Hernquist [35] sphere 

*W = -^^, 9 

r + Ohr 

with = 0.5 and Mhr = 0.2. In this case, the Plummer and Hernquist spheres 
play the role of a dark halo and bulge, respectively. The total relative mass of the 
spherical subsystem within four exponential disk parameters is approximately unity 
(Msph(4/i)/Mdisk(4/i) ~ 1). The duration of a single iteration is ti = 15; the thickness 
of the cylindrical layers into which we subdivided the model for the transfer of the ve- 
locity distribution function is = 0.1 (i.e., a less coarse dividing than in the first case), 
and the number of particles is = 5 ■ 10^. 

The models constructed using the iterative method have the following properties. 

1. The models prove to be close to equilibrium. They preserve both their struc- 
tural and dynamical parameters well during the initial stages of their evolution 
(Figs. [9] and [To]). It goes without saying that different values of pdisk and $ext yield 
models with different degrees of closeness to equilibrium. For example, the models 
of a thin stellar disk without an external potential (with the density distribution 
([21 and h = 3.5, Zq = 0.3, = 1, -Rmax = 14) are appreciably nonequilibrium. 
However, in any case, the resulting models are substantially closer to equilibrium 
than those constructed using the Hernquist technique (see Section [2]). We consider 
possible reasons for deviations of the iterative models from equilibrium below. 
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Figure 7: Comparison of the moments of the velocity distribution function for the iterative 
model of the first family with the solutions of the Jeans equations ([T]). In the left-hand 
plot, the solid curve (model) shows the mean azimuthal velocity for the iterative model, 
and the dashed curve (Jeans equation) shows the same quantity computed proceeding 
from the Jeans equation [the first equation of ([T])], where we adopted aR and from the 
iterative model. In the middle plot, the solid curve shows the azimuthal velocity dispersion 
for the iterative model, and the dashed curve the same quantity computed using the Jeans 
equations [the second equation of ([T])], where we adopted and aR from the iterative 
model. In the right-hand plot, the solid curve shows the dispersion of the stellar velocities 
in the vertical direction, and the dashed curve the same quantity computed using the 
Jeans equations [the third equation of ([T])]. We used a model of the family presented in 
Fig. O (marked by a cross). 




Figure 8: Same as Fig. [7] for the iterative model of the second family. We used a model 
of the family shown in Fig. [6] (marked by a cross) . 
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Figure 9: Initial stages of the evolution of an iterative model in the first family. A model 
of the family shown in Fig. [5] (marked by a cross) is used. The top snapshots show a face- 
on view for several times, with the shades of gray corresponding to the logarithm of the 
particle number per pixel. The plots in the central and bottom parts of the figure show 
the dependences of various disk parameters on the cylindrical radius R for various times: 
the number of particles in concentric cylindrical layers n, twice the median of 1^1, 2zi/2 
(which characterizes the disk thickness and is close to zq for the distribution ([2])), and 
the four moments of the velocity distribution function -u^, a^, a^, and az- N = 25 000, 
e = 0.02, and dt = 0.01. The disk thickness can be seen to have changed slightly. This is 
probably due to the assumption that the 2;-velocity distribution function is isothermal. 
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Figure 10: Same as Fig. [9] for an iterative model of the second family {N = 50000, 
e = 0.02, dt = 0.01). A model of the family shown in Fig. [6] (marked by a cross) is used. 
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2. The first four moments of the velocity distribution function fairly accurately obey 
the equilibrium Jeans equations ([T]) (Figs. Hand [H]). This means that, if we replaced 
the assumed dependence aji{R) in the Hernquist method (see Section [2]) with the 
dependence obtained from the iterations, the resulting models will be very similar 
to the iterative models, i.e., they essentially be in equilibrium. In other words, given 
a more correct assumption about the dependence aR^R), the Hernquist technique 
would allow the construction of close-to-equilibrium models. 

3. All models in some family of iterative models have the same vertical- velocity dis- 
persion profile az{R) (Figs. [5] and [6]). This means that the equilibrium condition in 
the vertical direction is completely decoupled from the disk equilibrium condition 
in the plane. This was expected, since, of all the moments of the velocity distribu- 
tion function, only cr^ appears in the Jeans equation describing equilibrium in the 
vertical direction [the last equation of ([I])]. 

4. The profile of the radial stellar-velocity dispersion, a^, differs strongly from the 
commonly adopted exponential profile (Figs. [5] and E]). 

4.2.3 Specific features of the iterative disk models. Further development of 
the technique 

Let us now turn to some specific features of our iterative models. 

1. Our stellar disk model is not a very good choice. The density profile in this 
model has a sharp cutoff. Naturally, such a structure cannot be in equilibrium, and 
the disk must have a smooth cutoff. We can see in Figs. [9] and [10 (the left-hand plot 
in the middle row) how this sharp cutoff is disrupted. We conclude that we must 
choose an initial stellar-disk model with a smooth density decrease at R values close 

to -Rmax- 

2. As is clear from Figs. [9] and [ID] (middle plot in the middle row), the thickness 
of the stellar disk varies somewhat during the initial stages of evolution. These 
variations can be seen in almost all the iterative models, and are more pronounced 
in models that are hotter in the disk plane and less pronounced in models that are 
cooler in the disk plane (here, hotter means that a higher fraction of the kinetic 
energy is contained in random motions). This nonequilibrium state of the iterative 
models is likely due to the assumptions about the velocity distribution function 
adopted when implementing our method, in particular our assumption that the 
velocity distribution function does not depend on z (the model is isothermal in 
the z direction). In addition, the assumption that vrv^ = may not be valid 
in the central regions. Future analyses should apply the iterative method without 
these assumptions, first and foremost, without the isothermal assumption, especially 
since this method is easy to generalize in this way. In the technique described 
here, we subdivided the model into cylindrical layers when transfering the velocity 
distribution function from the previous iteration step. Each of these cylindrical 
layers can be subdivided into layers along the z axis, and the moments of the 
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velocity distribution function can be computed separately in each layer (i.e., the 
model can be subdivided in both the R and z planes). Note, however, that such a 
modification of the iterative method requires a far greater number of particles than 
we have used so far. 

3. The iterative method cannot yield an equilibrium model that is unstable on time 
scales shorter than the time for one iteration. For example, the method cannot 
produce a very cold model. As is evident in Figs. [5] and El the family of iterative 
models is bounded on the side of cold models. It is not possible to construct a very 
cold model, even when is fixed - the iterations still fail to converge. 

4.2.4 Assumption of uniqueness of the family of models. New applications 
for iterative models 

Thus, the technique described above can be used to construct a one-parameter family 
of close-to-equilibrium models for a given pdisk and $ext- This family is parametrized 
by the fraction of kinetic energy contained in random motions (e.g., this parameter can 
have the form of the total angular momentum of the disk with respect to the z axis, L^). 
The question naturally arises in this case of whether there are other equilibrium models 
besides those in the family of iterative disk models? We can only postulate that there can 
be only one equilibrium disk model. Such a hypothesis can be formulated follows. There 
can be at most one equilibrium model (one equilibrium distribution function) for a given 
Pdisk{R, z), ^ext{R, z) with a fixed fraction of kinetic energy contained in random motions 
(e.g., for fixed L^). 

We can neither prove nor disprove this statement, and can only assume that it is 
true. This assumption expands the scope of application of the iterative models. We 
can compare iterative models to observational data in order to derive constraints on 
unobserved parameters of galaxies. For example, our statement that the dependence of 
cT/j on R for the iterative models is far from exponential requires additional explanations. 
More precisely, this dependence is far from exponential for the pdisk and $ext considered. 
It is not ruled out that an external potential can be selected so as to make the dependence 
of ctr on R exponential p] (see also Section [3l where we discuss the observational data 
for this dependence). When available, GAIA data on the velocity field in our Galaxy 
can be used to try to constrain the mass distribution in the dark halo of the Galaxy via 
iterative models (by selecting an external potential that makes the equilibrium stellar disk 
embedded in this potential have the observed velocity field). 

5 Conclusions 

We have proposed here a new iterative approach to constructing equilibrium A^-body 
models with a given density distribution. The main idea of this approach is the following. 
At the first stage, a model is constructed using some approximate method. The model is 
then allowed to adjust to the equilibrium state while its density distribution is kept fixed, 
as well as the required parameters of the velocity distribution, if necessary. 
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We used our iterative approach to construct two types of models. The first type have 
the form of isotropic and spherically symmetric systems. We used the iterative approach 
to construct such systems in order to test whether the approach indeed enables the con- 
struction of close-to-equilibrium models. The second type of models were axisymmetric 
stellar disks in an external potential. In both cases, the numerical models constructed 
using the iterative approach were very close to equilibrium. The iterative models for the 
stellar disks were much closer to equilibrium than models constructed based on the Jeans 
equations (see Section [2]). 

We also show that the hypothesis that the radial velocity dispersion is proportional 
to the vertical velocity dispersion in spiral galaxies {an oc az) may be incorrect. First, 
this assumption lacks clear observational evidence to support it. Second, it is not true 
for the iterative models of stellar disks that we have constructed. Finally, the assumption 
that cT/j oc (T2 is the main reason why the stellar disks constructed using the technique of 
Hernquist [H] are fairly far from equilibrium. To prove this last statement, recall that, since 
the iterative models of stellar disks obey the Jeans equations, if we replace in Hernquist's 
technique the relation a/j oc cr^ with the dependence ctr^R) = a}^{R), where (j}i{R) is the 
radial dispersion profile for an iterative model, the resulting models should be as close to 
equilibrium as the iterative models. 

The proposed method has a wide range of applications, and can be used to construct 
models with various geometries. For example, it can be used to construct spherically sym- 
metric models with anisotropic stellar motions, as well as self-consistent multicomponent 
models of disk galaxies. Unlike the kinetic approach based on the distribution function, 
our method provides direct input data for A^-body simulations. 

Furthermore, the proposed method has purely astrophysical applications. For exam- 
ple, bringing the velocity field obtained in equilibrium iterative models into agreement 
with the velocity field in our Galaxy can enable determination of the parameters of the 
external potential, and, consequently, the dark halo. 
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